overview

The analysis below is an exploratory analysis of the coelacanth data set from Nikaido et. al.. It should be noted that the replicates are actually technical replicates and analyzing them using this method greatly underestimates the variance and as a result is not a valid population level inference, or even a valid sample specific inference. The final section (groupings samples together) covers a slightly better case, though still technically not what you should do if you have the opportunity to generate biological replicates.

Regardless, here we perform analysis of the tissue specific samples from Nikaido et. al. as an illustration of the types of analyses that you can do with sleuth.

loading libraries

First, we load the sleuth package. Next, we load cowplot which has some nice formatting modifications of the standard ggplot2 figures.

library('sleuth')
library('cowplot')

naive analysis using technical replicates

Below we load the sample to condition table which gives and explanation of the data.

sample_to_condition <- read.table('SraRunTable.txt', header = TRUE,
  stringsAsFactors = FALSE, sep = '\t')
head(sample_to_condition)
##    BioSample_s Center_Name_s
## 1 SAMD00004034           NIG
## 2 SAMD00004036           NIG
## 3 SAMD00004033           NIG
## 4 SAMD00004032           NIG
## 5 SAMD00004037           NIG
## 6 SAMD00004035           NIG
##                                                                         Library_Name_s
## 1       RNA from kidney of Tanzanian coelacanth Latimeria chalumnae, paired-end, 200bp
## 2   RNA from pelvic fin of Tanzanian coelacanth Latimeria chalumnae, paired-end, 200bp
## 3 RNA from pectoral fin of Tanzanian coelacanth Latimeria chalumnae, paired-end, 200bp
## 4      RNA from pharynx of Tanzanian coelacanth Latimeria chalumnae, paired-end, 200bp
## 5         RNA from gill of Tanzanian coelacanth Latimeria chalumnae, paired-end, 200bp
## 6  RNA from tail muscle of Tanzanian coelacanth Latimeria chalumnae, paired-end, 200bp
##     LoadDate_s MBases_l MBytes_l ReleaseDate_s       run SRA_Sample_s
## 1 Aug 07, 2012    19301    12048  Aug 07, 2012 DRR002302    DRS001630
## 2 Aug 07, 2012    19378    11642  Aug 07, 2012 DRR002303    DRS001631
## 3 Aug 07, 2012    16107    10114  Aug 07, 2012 DRR002304    DRS001632
## 4 Aug 07, 2012    25944    16419  Aug 07, 2012 DRR002305    DRS001633
## 5 Aug 07, 2012    25028    16030  Aug 07, 2012 DRR002306    DRS001634
## 6 Aug 07, 2012    15901     9638  Aug 07, 2012 DRR002307    DRS001635
##   Sample_Name_s
## 1  SAMD00004034
## 2  SAMD00004036
## 3  SAMD00004033
## 4  SAMD00004032
## 5  SAMD00004037
## 6  SAMD00004035
##                                                     sample_comment_s
## 1       RNA from kidney of Tanzanian coelacanth Latimeria chalumnae.
## 2   RNA from pelvic fin of Tanzanian coelacanth Latimeria chalumnae.
## 3 RNA from pectoral fin of Tanzanian coelacanth Latimeria chalumnae.
## 4      RNA from pharynx of Tanzanian coelacanth Latimeria chalumnae.
## 5         RNA from gill of Tanzanian coelacanth Latimeria chalumnae.
## 6  RNA from tail muscle of Tanzanian coelacanth Latimeria chalumnae.
##   sample_name_s Assay_Type_s AssemblyName_s BioProject_s Consent_s
## 1     DRS001630      RNA-Seq <not provided>    PRJDB1984    public
## 2     DRS001631      RNA-Seq <not provided>    PRJDB1984    public
## 3     DRS001632      RNA-Seq <not provided>    PRJDB1984    public
## 4     DRS001633      RNA-Seq <not provided>    PRJDB1984    public
## 5     DRS001634      RNA-Seq <not provided>    PRJDB1984    public
## 6     DRS001635      RNA-Seq <not provided>    PRJDB1984    public
##   InsertSize_l LibraryLayout_s          Organism_s Platform_s SRA_Study_s
## 1            0          PAIRED Latimeria chalumnae   ILLUMINA   DRP000627
## 2            0          PAIRED Latimeria chalumnae   ILLUMINA   DRP000627
## 3            0          PAIRED Latimeria chalumnae   ILLUMINA   DRP000627
## 4            0          PAIRED Latimeria chalumnae   ILLUMINA   DRP000627
## 5            0          PAIRED Latimeria chalumnae   ILLUMINA   DRP000627
## 6            0          PAIRED Latimeria chalumnae   ILLUMINA   DRP000627
##   g1k_analysis_group_s g1k_pop_code_s       source_s
## 1       <not provided> <not provided> <not provided>
## 2       <not provided> <not provided> <not provided>
## 3       <not provided> <not provided> <not provided>
## 4       <not provided> <not provided> <not provided>
## 5       <not provided> <not provided> <not provided>
## 6       <not provided> <not provided> <not provided>

A lot of these fields are redundant and hard to read. Let’s clean up the formatting a little bit.

stc <- dplyr::select(sample_to_condition, sample = run,
  center = Center_Name_s, tissue = Library_Name_s)
stc <- dplyr::mutate(stc,
  tissue = sapply(strsplit(tissue, ' from '), '[[', 2))
stc <- dplyr::mutate(stc,
  tissue = sapply(strsplit(tissue, ' of'), '[[', 1))

Now, let’s add the path to the kallisto results:

stc <- dplyr::mutate(stc,
  path = file.path('results', sample, 'kallisto'))

Everything is a lot easier to read now:

head(stc)
##      sample center       tissue                       path
## 1 DRR002302    NIG       kidney results/DRR002302/kallisto
## 2 DRR002303    NIG   pelvic fin results/DRR002303/kallisto
## 3 DRR002304    NIG pectoral fin results/DRR002304/kallisto
## 4 DRR002305    NIG      pharynx results/DRR002305/kallisto
## 5 DRR002306    NIG         gill results/DRR002306/kallisto
## 6 DRR002307    NIG  tail muscle results/DRR002307/kallisto

Note that for every single sleuth run, there should be at least 2 columns of metadata:

You should also have additional columns that explain the data (experimental conditions, batch, etc.).

loading gene names

While not required, below we load the gene names using the biomaRt package. This will attach a gene name to each transcript and is also required for gene aggregation mode.

mart <- biomaRt::useMart(biomart = "ensembl",
  dataset = 'lchalumnae_gene_ensembl')

ttg <- biomaRt::getBM(
  attributes = c("ensembl_transcript_id", "ensembl_gene_id", "external_gene_name"),
  mart = mart)
ttg <- dplyr::rename(ttg,
  target_id = ensembl_transcript_id,
  ens_gene = ensembl_gene_id,
  ext_gene = external_gene_name)

If loading this table into sleuth, you must have a column titled target_id and another column that represents the gene (or other metadata).

Here is a sample of what it looks like:

head(ttg)
##            target_id           ens_gene ext_gene
## 1 ENSLACT00000000005 ENSLACG00000000004         
## 2 ENSLACT00000000150 ENSLACG00000000131         
## 3 ENSLACT00000000550 ENSLACG00000000487         
## 4 ENSLACT00000000654 ENSLACG00000000581         
## 5 ENSLACT00000000655 ENSLACG00000000582         
## 6 ENSLACT00000000656 ENSLACG00000000583

You might want to save this for later in case you don’t have internet access. Use readRDS() to load this file later.

saveRDS(ttg, file = 'ttg.rds')

prepping the object and loading data

We are now ready to load the data. Sleuth requires:

so <- sleuth_prep(stc, ~tissue + center, target_mapping = ttg,
  max_bootstrap = 30)
## reading in kallisto results
## ...................
## 
## normalizing est_counts
## 17842 targets passed the filter
## normalizing tpm
## merging in metadata
## normalizing bootstrap samples
## summarizing bootstraps

Let’s look at some of the diagnostics.

sleuth_live(so)

A few things you might want to look at are:

plot_pca(so, text_labels = TRUE, color_by = 'tissue')

Notice that tissue specific samples cluster together.

We can also change the variable is colored by:

plot_pca(so, text_labels = TRUE, color_by = 'center')

From this PCA it is pretty clear that there is a sequencing center specific effect.

You might also want to check out how the replicates did. Below are two of the gill samples sequenced at different centers:

plot_scatter(so, 'DRR002306', 'DRR002319')

Additionally, other things that you could look at that are interesting and useful:

There are plenty of other things to look at – please explore sleuth_live()!

fitting the models

Next, we fit the full model including a regressive for the batch (center).

so <- sleuth_fit(so)
## fitting measurement error models
## shrinkage estimation
## Adding missing grouping variables: `x_group`
## computing variance of betas

Next, we regress only against the batch (center) and call the model ‘center’. This allows us to condition on the center when testing.

so <- sleuth_fit(so, ~center, 'center')
## fitting measurement error models
## shrinkage estimation
## Adding missing grouping variables: `x_group`
## computing variance of betas
so <- sleuth_fit(so, ~1, 'intercept')
## fitting measurement error models
## shrinkage estimation
## Adding missing grouping variables: `x_group`
## computing variance of betas

For fun, we can also see what we get out when we ignore the sequencing center:

so <- sleuth_fit(so, ~tissue, 'tissue')
## fitting measurement error models
## shrinkage estimation
## Adding missing grouping variables: `x_group`
## computing variance of betas

Now, we can actually compared the models conditioning on center and center and tissue to find effects specific to tissue while conditioning on the center.

so <- sleuth_lrt(so, 'center', 'full')

Let’s take a look at things that are significant by the batch and sequencing center alone:

so <- sleuth_lrt(so, 'intercept', 'center')
so <- sleuth_lrt(so, 'intercept', 'tissue')

Now, we can interactively look at things some more with the tests included:

sleuth_live(so)

If we want to look at the results in R, we can load them using sleuth_results.

intercept_center <- sleuth_results(so, 'intercept:center', 'lrt')
intercept_tissue <- sleuth_results(so, 'intercept:tissue', 'lrt')
center_full <- sleuth_results(so, 'center:full', 'lrt')

Let’s take a look at the ‘significant’ things that come out:

level <- 0.05

nrow(dplyr::filter(intercept_center, qval < level))
## [1] 1528
nrow(dplyr::filter(intercept_tissue, qval < level))
## [1] 12882
nrow(dplyr::filter(center_full, qval < level))
## [1] 15410

We can remove the hits that are significant from the batch.

remove_batch <- dplyr::anti_join(
  dplyr::filter(center_full, qval < level),
  dplyr::filter(intercept_center, qval < level),
  by = 'target_id'
  )
remove_batch <- dplyr::arrange(remove_batch, qval, desc(test_stat))

Here, we take a look at some of the potential false positives if we had not removed them from batch.

possible_false_positives <- dplyr::semi_join(
  dplyr::filter(center_full, qval < level),
  dplyr::filter(intercept_center, qval < level),
  by = 'target_id'
  )
possible_false_positives <- dplyr::arrange(possible_false_positives, qval, desc(test_stat))

Let’s look at the top hits from both of these tests:

remove_batch %>% head
##            target_id test_stat         pval         qval       rss
## 1 ENSLACT00000004099  115.8770 2.342322e-23 4.177062e-19 122.68544
## 2 ENSLACT00000015097  112.7313 1.084159e-22 9.666905e-19  76.94338
## 3 ENSLACT00000016039  109.6389 4.884444e-22 2.316432e-18  66.61245
## 4 ENSLACT00000007913  109.5119 5.195832e-22 2.316432e-18  60.03172
## 5 ENSLACT00000004505  108.5717 8.209091e-22 2.510220e-18  51.59473
## 6 ENSLACT00000020776  108.4743 8.607201e-22 2.510220e-18  61.92421
##   sigma_sq     tech_var  mean_obs  var_obs sigma_sq_pmax smooth_sigma_sq
## 1 7.206998 9.792204e-03  8.386656 6.815894      7.206998       0.1521885
## 2 4.525990 9.124694e-05 11.445480 4.294710      4.525990       0.5304479
## 3 3.918318 6.133232e-05 11.281432 3.751509      3.918318       0.4894698
## 4 3.531157 1.209962e-04 10.703871 3.344419      3.531157       0.3722032
## 5 3.034674 3.098772e-04  9.523322 2.924275      3.034674       0.2249321
## 6 3.642549 5.180990e-05 11.648330 3.531637      3.642549       0.5867137
##   final_sigma_sq degrees_free           ens_gene ext_gene
## 1       7.206998            5 ENSLACG00000003616    TNNT1
## 2       4.525990            5 ENSLACG00000013194   mybpha
## 3       3.918318            5 ENSLACG00000014026    ttn.2
## 4       3.531157            5 ENSLACG00000006948         
## 5       3.034674            5 ENSLACG00000003974    MYH15
## 6       3.642549            5 ENSLACG00000018133    MYH7B
possible_false_positives %>% head
##            target_id test_stat         pval         qval      rss
## 1 ENSLACT00000012180  62.88042 3.081999e-12 4.942561e-11 3.673126
## 2 ENSLACT00000021477  56.26037 7.182632e-11 7.995498e-10 4.032439
## 3 ENSLACT00000007217  53.77458 2.331736e-10 2.242819e-09 3.164205
## 4 ENSLACT00000005300  52.62668 4.012402e-10 3.619280e-09 2.761997
## 5 ENSLACT00000006114  52.15186 5.021416e-10 4.406836e-09 5.432517
## 6 ENSLACT00000008461  51.87358 5.726648e-10 4.957442e-09 9.818850
##    sigma_sq     tech_var mean_obs   var_obs sigma_sq_pmax smooth_sigma_sq
## 1 0.2154044 0.0006618229 8.117402 0.3204511     0.2154044       0.1410221
## 2 0.2335699 0.0036323805 7.492555 0.3658917     0.2335699       0.1213775
## 3 0.1840495 0.0020802589 6.731672 0.3155036     0.1840495       0.1073984
## 4 0.1614772 0.0009931938 7.249908 0.2392036     0.1614772       0.1157234
## 5 0.3178490 0.0017108225 6.420179 0.5268688     0.3178490       0.1068910
## 6 0.5734458 0.0041335849 5.845879 1.0809316     0.5734458       0.1184288
##   final_sigma_sq degrees_free           ens_gene ext_gene
## 1      0.2154044            5 ENSLACG00000010642   PITPNB
## 2      0.2335699            5 ENSLACG00000018747    PSMB7
## 3      0.1840495            5 ENSLACG00000006350    PAIP2
## 4      0.1614772            5 ENSLACG00000004667    SUGT1
## 5      0.3178490            5 ENSLACG00000005382    S100B
## 6      0.5734458            5 ENSLACG00000007433  CEP170B

The first few hits look very tissue specific:

plot_bootstrap(so, remove_batch$target_id[1], color_by = 'tissue')

plot_bootstrap(so, remove_batch$target_id[5], color_by = 'tissue')

plot_bootstrap(so, remove_batch$target_id[10], color_by = 'tissue')

We can also look at possible false positives that are possibly associated with batch as well:

plot_bootstrap(so, possible_false_positives$target_id[1], color_by = 'center')

plot_bootstrap(so, possible_false_positives$target_id[5], color_by = 'center')

plot_bootstrap(so, possible_false_positives$target_id[10], color_by = 'center')

You can more easily interactively look at these results under the sleuth live tab analyses -> transcript view

sleuth_live(so)

Let’s see if we can summarize the top hits with some figures. First, let’s pull out the top hits.

top_100 <- dplyr::filter(center_full, qval < level) %>%
  head(100)

Next, we have to convert the observations to a matrix object for most heatmap functions:

tmp <- sleuth_to_matrix(so, 'obs_norm', 'est_counts')$data
tmp <- tmp[top_100$target_id, ]

Next, let’s rename the columns so that we can easily identify the sample types.

colnames(tmp)
##  [1] "DRR002307" "DRR002311" "DRR002313" "DRR002310" "DRR002304"
##  [6] "DRR002314" "DRR002320" "DRR002303" "DRR002309" "DRR002316"
## [11] "DRR002306" "DRR002315" "DRR002318" "DRR002308" "DRR002302"
## [16] "DRR002317" "DRR002319" "DRR002312" "DRR002305"
stc <- stc[match(colnames(tmp), stc$sample), ]
colnames(tmp) <- paste0(stc$tissue, '_', stc$center)

Finally, we can make the heatmap.

heatmap(log(tmp + 1))

groupings samples together

Notice, the top differentially expressed things cluster by (pectoral fin, pelvic fin) and (pharynx, tail muscle). This makes sense because the first two samples are fin types and the second two samples are muscle types.

While we don’t exactly have biological replicates, we can model some of the variance within similar tissues by pooling them together.

NOTE: This will sort of tell us about tissue specific expression within this fish, but not at the population level. In short, this is basically an exploratory analysis.

replicate_mapping <- data.frame(
  tissue = c('pectoral fin', 'pelvic fin', 'pharynx', 'tail muscle'),
  grouping = c('fin', 'fin', 'muscle', 'muscle'),
  stringsAsFactors = FALSE)

stc_grouped <- dplyr::inner_join(stc, replicate_mapping, by = 'tissue')

For the simple analysis, let’s remove the batch effect by simply ignoring everything from that center.

to_keep <- stc_grouped %>% dplyr::filter(center == 'NIG')

Also, let’s drop the technical replicates the we are modeling variance between different cells and not just from the same library prep but a different run.

set.seed(43)
drop_replicates <- to_keep %>%
  dplyr::group_by(tissue) %>%
  dplyr::sample_n(1)

Let’s run sleuth_prep on this new grouping.

so_pool <- sleuth_prep(drop_replicates, ~grouping, target_mapping = ttg,
  max_bootstrap = 30)
## reading in kallisto results
## ....
## 
## normalizing est_counts
## 19081 targets passed the filter
## normalizing tpm
## merging in metadata
## normalizing bootstrap samples
## summarizing bootstraps

Let’s perform tests now against the full model (grouped tissues) versus the intercept only model.

so_pool <- sleuth_fit(so_pool)
## fitting measurement error models
## shrinkage estimation
## Adding missing grouping variables: `x_group`
## computing variance of betas
so_pool <- sleuth_fit(so_pool, ~1, 'intercept')
## fitting measurement error models
## shrinkage estimation
## Adding missing grouping variables: `x_group`
## computing variance of betas
so_pool <- sleuth_lrt(so_pool, 'intercept', 'full')

pool_result <- sleuth_results(so_pool, 'intercept:full', 'lrt')
pool_result <- dplyr::arrange(pool_result, qval, desc(test_stat))

Let’s look at a heatmap of the top 100 things again.

pool_top_100 <- pool_result[1:100, ]

tmp <- sleuth_to_matrix(so_pool, 'obs_norm', 'est_counts')$data
# tmp <- sleuth_to_matrix(so_pool, 'obs_norm', 'tpm')$data
tmp <- tmp[pool_top_100$target_id, ]
drop_replicates <- drop_replicates[match(colnames(tmp), drop_replicates$sample), ]

colnames(tmp) <- paste0(drop_replicates$grouping)

heatmap(log(tmp + 1))

Let’s look at the top 15 results.

pool_result %>% head(15)
##             target_id test_stat         pval        qval      rss
## 1  ENSLACT00000018951  24.20607 8.655914e-07 0.002167192 14.15779
## 2  ENSLACT00000001320  23.60884 1.180443e-06 0.002167192 15.69639
## 3  ENSLACT00000010421  23.33370 1.361924e-06 0.002167192 13.89557
## 4  ENSLACT00000018864  23.27800 1.401944e-06 0.002167192 19.58818
## 5  ENSLACT00000005548  23.15919 1.491288e-06 0.002167192 11.64981
## 6  ENSLACT00000007913  23.10647 1.532742e-06 0.002167192 18.39947
## 7  ENSLACT00000004620  23.02197 1.601606e-06 0.002167192 15.29150
## 8  ENSLACT00000016667  23.01628 1.606350e-06 0.002167192 10.25620
## 9  ENSLACT00000005838  22.94906 1.663514e-06 0.002167192 11.64983
## 10 ENSLACT00000015097  22.31428 2.314789e-06 0.002167192 21.81623
## 11 ENSLACT00000018278  22.27213 2.366161e-06 0.002167192 10.83293
## 12 ENSLACT00000001981  22.25980 2.381396e-06 0.002167192 80.08362
## 13 ENSLACT00000000023  22.20117 2.455235e-06 0.002167192 14.54017
## 14 ENSLACT00000020580  22.12930 2.548897e-06 0.002167192 15.15890
## 15 ENSLACT00000005635  22.04088 2.669053e-06 0.002167192 12.54392
##     sigma_sq     tech_var  mean_obs   var_obs sigma_sq_pmax
## 1   4.718848 4.157621e-04  8.686350  4.719264      4.718848
## 2   5.231696 4.333184e-04  8.974983  5.232129      5.231696
## 3   4.630171 1.685954e-03  7.317912  4.631857      4.630171
## 4   6.527491 1.901381e-03  7.098046  6.529392      6.527491
## 5   3.883071 1.994719e-04  9.409498  3.883271      3.883071
## 6   6.133040 1.166769e-04 11.017716  6.133157      6.133040
## 7   5.097096 7.070713e-05 10.474561  5.097166      5.097096
## 8   3.418459 2.745971e-04  9.051549  3.418734      3.418459
## 9   3.882924 3.533803e-04  9.075228  3.883277      3.882924
## 10  7.272007 6.876526e-05 11.575784  7.272076      7.272007
## 11  3.610033 9.437372e-04  7.749520  3.610976      3.610033
## 12 26.693863 6.778831e-04  3.779807 26.694541     26.693863
## 13  4.846681 4.315931e-05 11.173876  4.846724      4.846681
## 14  5.049001 3.966470e-03  6.485851  5.052967      5.049001
## 15  4.178772 2.535907e-03  6.889045  4.181308      4.178772
##    smooth_sigma_sq final_sigma_sq degrees_free           ens_gene
## 1        0.1784171       4.718848            1 ENSLACG00000016564
## 2        0.1907583       5.231696            1 ENSLACG00000001174
## 3        0.1432852       4.630171            1 ENSLACG00000009108
## 4        0.1423120       6.527491            1 ENSLACG00000016490
## 5        0.2135253       3.883071            1 ENSLACG00000004890
## 6        0.3632446       6.133040            1 ENSLACG00000006948
## 7        0.2980085       5.097096            1 ENSLACG00000004079
## 8        0.1943840       3.418459            1 ENSLACG00000014586
## 9        0.1955370       3.882924            1 ENSLACG00000005138
## 10       0.4527856       7.272007            1 ENSLACG00000013194
## 11       0.1503622       3.610033            1 ENSLACG00000015985
## 12       0.5241814      26.693863            1 ENSLACG00000001757
## 13       0.3857136       4.846681            1 ENSLACG00000000018
## 14       0.1527309       5.049001            1 ENSLACG00000017964
## 15       0.1433523       4.178772            1 ENSLACG00000004966
##           ext_gene
## 1                 
## 2          ST8SIA4
## 3                 
## 4  pgp (1 of many)
## 5           FBXO32
## 6                 
## 7             LDB3
## 8                 
## 9          SLC38A2
## 10          mybpha
## 11         APOBEC2
## 12           PSME4
## 13           MYOM2
## 14            FSD2
## 15          CACNG6

We notice, that the gene MYH15 is in this list. This is a Myosin gene, associated with slow twitch muscle tissue in other species.

plot_bootstrap(so_pool, 'ENSLACT00000004505', color_by = 'grouping')

Like in our other analysis, you can also explore the results using sleuth_live():

sleuth_live(so_pool)

session info

sessionInfo()
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] shiny_0.14.1  cowplot_0.6.3 sleuth_0.28.1 dplyr_0.5.0   ggplot2_2.1.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.7          formatR_1.4          plyr_1.8.4          
##  [4] bitops_1.0-6         tools_3.3.1          zlibbioc_1.18.0     
##  [7] biomaRt_2.28.0       digest_0.6.10        evaluate_0.10       
## [10] jsonlite_1.1         RSQLite_1.0.0        tibble_1.2          
## [13] gtable_0.2.0         rhdf5_2.16.0         DBI_0.5-1           
## [16] yaml_2.1.13          parallel_3.3.1       knitr_1.14          
## [19] stringr_1.1.0        S4Vectors_0.10.3     IRanges_2.6.1       
## [22] stats4_3.3.1         grid_3.3.1           data.table_1.9.6    
## [25] Biobase_2.32.0       R6_2.2.0             AnnotationDbi_1.34.4
## [28] XML_3.98-1.4         rmarkdown_1.1        tidyr_0.6.0         
## [31] reshape2_1.4.2       magrittr_1.5         scales_0.4.0        
## [34] htmltools_0.3.5      BiocGenerics_0.18.0  assertthat_0.1      
## [37] xtable_1.8-2         mime_0.5             colorspace_1.2-6    
## [40] httpuv_1.3.3         labeling_0.3         stringi_1.1.2       
## [43] RCurl_1.95-4.8       lazyeval_0.2.0       munsell_0.4.3       
## [46] chron_2.3-47